function ivol = blsImpVolPut( S, K, r, T, price, q, tolerance )

xn = 0.4*ones(length(S),1);

for i = 1:50
    temp1 = ( log(S./K)+r-q )./( xn.*sqrt(T) );
    d1 = temp1 + xn.*sqrt(T)/2;
    d2 = temp1 - xn.*sqrt(T)/2;

    phi1 = exp( -d1.^2 / 2 ) / sqrt(2*pi);

    temp2 = xn - 1./sqrt(T).*( (exp(-r.*T).*K.*normcdf(-d2) - price)./(S.*exp(-q*T).*phi1) - normcdf(-d1)./phi1 );
    
    if max(abs(temp2-xn)) < tolerance
        ivol = temp2;
        return;
    else
        xn = temp2;
    end
end

ivol = NaN*ones(length(S),1);